Merger of black hole-neutron star binaries: nonspinning black hole case 
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We perform a simulation for merger of a black hole (BH)-neutron star (NS) binary in full general 
relativity preparing a quasicircular state as initial condition. The BH is modeled by a moving 
puncture with no spin and the NS by the F-law equation of state with F = 2. Corotating velocity 
field is assumed for the NS. The mass of the BH and the rest-mass of the NS are chosen to be 
w 3.2M and » 1.4M Q with relatively large radius of the NS » 14 km. The NS is tidally disrupted 
near the innermost stable orbit but ~ 80% of the material is swallowed into the BH with small disk 
mass ~ 0.3Mq even for such small BH mass ~ 3M . The result indicates that the system of a BH 
and a massive disk of ~ Mq is not formed from nonspinning BH-NS binaries, although a disk of 
mass ~ Q.IMq is a possible outcome. 
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Merger of black hole (BH)-neutron star (NS) binaries is 
one of likely sources of kilo- meter size laserinterferometric 
gravitational wave detectors. Although such system has 
not been observed yet, statistical studies based on the 
stellar evolution synthesis suggest that the merger will 
happen more than 10% as frequently as the merger of 
binary NSs [l], Thus, the detection of such system 
will be achieved by laserinterferometers in near future. 

According to a study based on the tidal approxima- 
tion (which is referred to as a study for configuration 
of a Newtonian star in circular orbits around a BH in 
its relativistic tidal field; e.g., [!, 0, H, [1, the fate 
is classified into two cases, depending on the mass ratio 
q = Mns/-^bH) where Mbh and Mns denote the masses 
of BH and NS, respectively. For q <, q c , the NS of radius 
R will be swallowed into the BH horizon without tidal 
disruption before the orbit reaches the innermost stable 
circular orbit (ISCO) @, H, while for q ^ q c , NS may 
be tidally disrupted before plunging into BH. Here, the 
critical value of q c depends on the BH spin and equation 
of state (EOS) of NS, and for the nonspinning case with 
stiff EOSs, q c re 0.3-0.35(i?/5M NS )- 3/2 @. (Throughout 
this paper, we adopt the geometrical units c = G = 1.) 

The second case has been studied with great interest 
because of the following reasons, (i) Gravitational waves 
at tidal disruption will bring information about the NS 
radius since the tidal disruption limit depends sensitively 
on it |8|. The relation between the mass and the radius 
of NSs may be used for determining the EOS of high 
density matter 0. (ii) Tidally disrupted NSs may form 
a massive disk of mass ~ 0.1— 1M© around the BH if 
the tidal disruption occurs outside the ISCO. Systems 
consisting of a BH and a massive, hot disk have been 
proposed as one of likely sources for the central engine 
of gamma-ray bursts (GRBs) with a short duration [lCj , 
and hence, merger of low-mass BH and NS is a candidate. 

However, the scenario based on the tidal approxima- 
tion studies may be incorrect since gravitational radia- 
tion reaction and gravitational effects of NS to the orbital 
motion are ignored. Radiation reaction shortens the time 
available for tidally disrupting NSs. The gravity of NS 



could increase the orbital radius of the ISCO and hence 
the critical value of the tidal disruption, q c , may be larger 
in reality. Miller ll| estimates these ignored effects and 
suggests that NSs of canonical mass and radius will be 
swallowed into BH without tidal disruption. Moreover, 
NSs are described by the Newtonian gravity in the tidal 
approximation. If we treat it in general relativity, the 
gravity is stronger and hence tidal disruption is less likely. 

Tidal disruption of NSs by a BH has been investigated 
in the Newtonian [13 ] and ap proximately general rela- 
tivistic (GR) simulation [3, [Hi- However, a simulation in 
full general relativity is required (see [14[ for an effort). 
In this paper, we present our first results for fully GR 
simulation, performed by our new code which has been 
improved from previous one (l5l flfl ] ; we enable to han- 
dle orbiting BHs adopting the moving puncture method 
recently developed [ItJ (see also [l8| for detailed calibra- 
tion). As the initial condition, we prepare a quasicircular 
state computed in a new formalism described below. In 
this paper, we focus on whether NSs of realistic mass 
and radius is tidally disrupted to form a massive disk 
around nonspinning BHs. We will illustrate that a disk 
with mass ~ Mq is an unlikely outcome for plausible 
values of NS mass and radius although a disk of mass of 
O(O.1M ) is possible. 
Formalism for a quasicircular state 

Three groups have worked in computing quasicircular 
states of BH-NS binaries [TH, HtJ Hl[ . Here, we propose a 
new method for computing accurate quasicircular states 
that can be used for numerical simulation in the moving 
puncture framework (l7| . 

Even just before the merger, it is acceptable to assume 
that BH-NS binaries are in a quasicircular orbit since 
the time scale of gravitational radiation reaction is a few 
times longer than the orbital period. Thus, we assume 
the presence of a helical Killing vector around the mass 
center of the system, £^ = (<9t) M 
orbital angular velocity f2 is constant 

Next, we assume that the NS is corotating around the 
mass center of the system. Irrotational velocity field is 
believed to be more realistic for BH-NS binaries [22| but 
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we do not choose it in this paper. The NSs in corotation 
are more subject to tidal disruption than in irrotation 
(e.g., [![) since the outer part has larger angular momen- 
tum than the inner part and thus outward mass ejection 
more easily occurs. If a corotating NS is stable against 
tidal disruption, this will be also the case for the irrota- 
tional NS of the same mass and radius. 

The assumption of corotating velocity field in the he- 
lical symmetric spacetime yields the first integral of the 



Euler equation, h 



= const, where h is specific en- 



thalpy defined by l + e + P/p, and e, P, and p are specific 
internal energy, pressure, and rest-mass density, respec- 
tively. In the present work, we adopt the T-law EOS 
with r = 2; P = pe = up 2 with n an adiabatic constant. 
u^ 1 denotes the four velocity and u l its time component. 
Assumption of corotation implies u M = u 1 ^ . 

For a solution of geometric variables of quasicircular 
orbits, we adopt the conformal flatness formalism for 
three-geometry. In this formalism, the solution is ob- 
tained by solving Hamiltonian and momentum constraint 
equations, and an equation for the time slicing condition 
which is derived from K k k = where Kij is the extrinsic 
curvature and K k k its trace (23[. Using the conformal 
factor %l> and the rescaled tracefree extrinsic curvature 
A? = ip^K^ , these equations are respectively written 

= -2npg$* - ll/AM- 7 , (1) 
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where A denotes the flat Laplacian, pn = ph(au 1 ) 2 — P, 
Ji = phm and S = ph^au 1 ) 2 — 1] + 3P. Here, a is the 
lapse function and $ is defined by $ = aip. 

We solve these equations in the framework of the punc- 
ture BH [It], H3, • Assuming that the puncture is lo- 
cated at rp, we set ip and $ 

' ' V, 
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where Mp and C are positive constants, and tbh 
(x BH = x k — Xp). Then elliptic equations for <p and rj are 
derived. The constant Mp is arbitrarily given, while C 
is determined from the virial relation (e.g., [26[) 

I diQdS* = - I d l ipdS i = 2ttM, (5) 

J r — >oo Jr — >oo 

where M is the ADM mass. The mass center is deter- 
mined from the condition that the dipole part of ip at 
spatial infinity is zero. 

Equation @ is rewritten setting 

i y (= A+Sfr) = W id + Wj,i - pijS^Wkj + Kf v (6) 

where Kfj denotes the weighted extrinsic curvature asso- 
ciated with linear momentum of a puncture BH; 

K ij = o^T~ \ n i p j + n 3 P i + ( n i n 3 ~ Sij)P k nkj- (7) 
BH V / 



Here, n k = n k = ^bhAbh- Pi denotes linear momentum 
of the BH, determined from the condition that the total 
linear momentum of system should be zero; 
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The RHS of Eq. © denotes the total linear momentum 
of the companion NS. Then, the total angular momentum 
of the system is derived from 



J = 



J J v il> 6 d 3 x + e zjk 46 kl Pi. 



The elliptic equation for Wi{= W l ) is 



(9) 
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Denoting W, = 7P>i — (x,i + B k ^x k ) where \ an d Pi are 
auxiliary functions [27|, Eq. (fTT))) is decomposed as 



AB, t = n,J^ b , A X = -TTJ iX l i> 



ijfi 



(11) 



Computing BH-NS binaries in a quasicircular orbit re- 
quires to determine the shift vector even in the puncture 
framework. This is because has to be obtained [it is 
derived from u k = 6 k iU t ip 4 '(v' 1 + fJ 1 ) where v l = £lip 1 }. The 
relation between (3 l and Aij is written 



5 jk dif3 k + 5 lk djf3 k 
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Operating S^di, an elliptic equation is derived 

A/3 1 + Istodkdjp* = 2d j {atp- 6 )A ij + 16naJj5 ij , (13) 
o 

which is solved in the same manner as that for Wi. 

We have computed several models of quasicircular 
states and found that the relation between Q and J ap- 
proximately agrees with the 3rd Post-Newtonian relation 
|28|. This makes us confirm that this approach is a fair 
way for preparing quasicircular states. We also found 
that in this method, the shift vector at r = rp automat- 
ically satisfies the condition j3 v = — f2 within the error of 
a few %. This implies that the puncture is approximately 
guaranteed to be in a corotating orbit in the solution. 
Numerical results for quasicircular states 

BH-NS binaries in quasicircular orbits have been com- 
puted for a wide variety of models with q = M*/Mbh ~ 
0.3-0.5 where M* denotes baryon rest-mass of the NS. In 
the present work, the compactness of spherical NSs with 
rest-mass M* is chosen to be « 0.14. In the T-law EOSs, 
the mass and radius of NS are rescaled by changing the 
value of k: In the following we fix the unit by setting that 
M* = IAMq. In this case, R « 13.8 km, the gravita- 
tional mass is 1.3OM , and M^/n 1 / 2 = 0.147. Accord- 
ing to theories for NSs based on realistic nuclear EOSs 
Hg, the radius of NS of M NS w 1.4M is 11-13 km. 
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FIG. 1: Snapshots of the density contour curves for p in the equatorial plane for model A. The solid contour curves are drawn for 
p = (l + 2i) x 10 14 g/cm 3 (i = 1, 2, 3) and for lO 14 " - 51 g/cm 3 (i = 1 ~ 8). The maximum density at t = is « 7.2 x 10 14 g/cm 3 . 
The blue, cyan, magenta, and green curves denote 10 14 , 10 12 , 10 11 , and 10 10 g/cm 3 , respectively. Vectors indicate the velocity 
field (v x ,v y ), and the scale is shown in the upper right-hand corner. The thick (red) circles are apparent horizons. 



TABLE I: Parameters of quasicircular states. Mass parameter 
of puncture, mass of BH, rest-mass of NS, mass and radius 
of NS in isolation, total mass of the system, nondimension 
angular momentum parameter, orbital period in units of M, 
and compactness of the system defined by C = (Aff2) 2//3 . 
Mass of BH is computed from the area of the apparent horizon 
A as (A/167r) 1//2 . Mass is shown in units of Mq. 
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A 3.13 


3.21 


1.40 


1.30 


13.8 


4.47 


0.729 


119 
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B 3.13 


3.21 


1.40 


1.30 


13.8 


4.47 


0.720 


110 


0.150 



Thus the radius chosen here is slightly larger than that 
of realistic NSs and is more subject to tidal disruption. 

In computation, we focus only on the orbit of slightly 
outside of ISCO. Table I shows the quantities for selected 
quasicircular states with q w 0.4; model A is that used for 
the following numerical simulation and model B is very 
close to the tidal dispruption limit of approximately the 
same mass as that of model A, showing that the model 
A has an orbit slightly outside the tidal disruption limit. 

The tidal approximation studies suggest that for q <; 
q« = 0.35(i?/5M NS )- 3 / 2 [(M BH ^)" 1 /6 3/2 ], NSs with T = 
2 could be tidally disrupted by a nonspinning BH 0. 
Here, ^ = M^/6 3/2 is the angular velocity of the ISCO 
around nonspinning BHs. For model B, q* m 0.32, 
and hence, q > <7*. According to the tidal approxima- 
tion studies 0, H|, such NS should be unstable against 
tidal disruption. Nevertheless, such equilibrium exists, 
proving that the tidal disruption limit in the frame- 
work of the tidal approximation does not give correct 
answer. Our studies indicate that the critical value g* is 
« 0.43(i?/5M NS r 3/2 [(M B Htt)~76 3/2 ]; tidal disruption 
of NS is much less likely than in the prediction by the 
tidal approximation Q . For the typical NS of radius 



will be necessary for (A/bh^)~ 2 ^ 3 > 6; this implies that 
canonical NSs will not be tidally disrupted outside ISCO 
by most of nonspinning BHs of mass larger than ~ 3M . 
Tidal disruption occurs only for NSs of relatively large 
radius and only for orbits very close to ISCO. 
Simulation for merger 

Even if tidal disruption of an NS occurs near ISCO, 
a massive disk may be formed around BH. To investi- 
gate this possibility we perform a numerical simulation 
adopting model A. 

For the simulation, we initially reset the lapse (i.e., 
<I>) since the relation a > should hold. In the present 
approach we give $ at t — by the following equation 



$ = 



n - 



1 + 0.1X 4 



1 + El=i^ ro + L1 ^ 



(14) 



where X = C/Vbh- Then, a = only at puncture and 
otherwise a > 0. Furthermore, for tbh > C, the values 
of <I> quickly approch to those of the quasicircular states. 

The numerical code for hydrodynamics is the same 
as that for performing merger of NS-NS binaries (high- 
resolution central scheme) 16]. However, we change 
equations for a, f3 l , and ip, and numerical scheme of 
handling the transport terms of evolution equations for 
geometries. For a and (3 l we solve 



{d t -^di)\na = 
d t (? = 0.75 f^F 3 



- AtdtFj) 



(15) 
(16) 



where 7^ is the conformal three-metric and Fi — 
5 s djjik- At denotes the time step for the simulation 



R ~ 5A/ns an d mass M\ 



NS 



1.4M Q , M BH £ 3.3M, 



and the second term in the RHS of Eq. (|16|) is intro- 
duced for stabilization. The equation for the conformal 
factor is also changed to 

d t tp~ 6 - d^- 6 f3 l ) = (aK k k - 2d l (3 l )tp- 6 , (17) 

since t(i diverges at the puncture [13] • 
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In addition, we have improved numerical scheme for 
the transport term of geometric variables (d t — j3 l di)Q 
where Q is one of the geometric variables: First, we 
rewrite this term to dtQ — di(Qj3 l ) + Qdi(3 l and then ap- 
ply the same scheme as in computing the transport term 
of the hydrodynamic equations to the second term (3rd- 
order interpolation scheme [15]). We have found that 
for evolving punctures, such high-resolution scheme for 
the transport term in the geometric variables are crucial. 
This is probably because of the fact that near punctures, 
some of geoemetric variables steeply vary and so is the 
term j3 l diQ. For other terms in the Einstein's equation, 
we use the 2nd-ordcr finite differencing as in 
(Note that in the case of nonuniform grid, 4-point finite 
differencing is adopted for Q^u since 3-point one is 1st- 
order.) After we performed most of runs, we repeated 
some of computations with a higher-order scheme as used 
in [13, EH]. With such scheme, convergent results are ob- 
tained with a relatively large grid spacing. However, the 
results are qualitatively unchanged and the extrapolated 
results (which are obtained in the limit of zero grid spac- 
ing) are approximately identical. 

We adopt the cell-centered Cartesian, (x,y,z), grid 
to avoid the situation that the location of punctures 
(which always stay in the z — plane) coincides with 
the grid location. The equatorial plane symmetry is as- 
sumed: the grid size is (27V, 27V, N) for x-y-z. Follow- 
ing [30|, we adopt a nonuniform grid; in the present 
approach, a domain of (2Nq,2Nq, Nq) grid zone is cov- 
ered with a uniform grid of the spacing Ax and out- 
side the domain, the grid spacing is increased accord- 
ing to £tanh[(i — No)/Ai]Ax where i denotes the i-th 
grid point in each direction. Nq, Ai, and £ are con- 
stants. Ax/Mp is chosen to be 1/8, 9/80, 1/10, 7/80, 
and 3/40. As shown in [l7l ]. such grid spacing can re- 
solve moving punctures. For (N, Nq, Ai, £, Ax/M p ), we 
choose (i) (160,105,30,4.5,1/8), (ii) (200,105,30,4.5,1/8), 
(iii) (200,105,30,5,1/8), (iv) (220,125,30,5,9/80), (v) 
(220,125,30,6,1/10), (vi) (220,140,30,7,7/80), and (vii) 
(220,150,9,3/40). For Ax = M P /8 and N = 160, we 
chose other values of Nq and Ai, and found that results 
depend weakly on them. The grid covers a cube of edge 
length 2L: For (i)-(vii), L/X = 0.46, 0.78, 0.83, 0.78, 
0.78, 0.65, and 0.59, respectively, where A is the wave- 
length of gravitational waves at t = 0. 

For a test, we performed a simulation for merger of two 
nonspinning BHs adopting the same initial condition as 
used in [llj • We focused particularly on the merger time 
and found that it varies with improving grid resolution. 
By the extrapolation, a true merger time is estimated. It 
is found that our result is ~ 19M and agrees with those 
of [TtJ (see 3l[ for our results). This indicates that our 



code can follow moving punctures as in [171 ]. 

Figure 1 shows evolution of contour curves for p and 
velocity vectors for v l in the equatorial plane together 
with the location of apparent horizons at selected time 
slices for (vii). Due to gravitational radiation reaction, 
the orbital radius decreases and then the NS is elongated 



(2nd panel). Because of the elongation, the quadrupole 
moment of the NS is amplified and the attractive force 
between two objects is strengthen [32j]. This effect accel- 
erates an inward motion and, consequently, the NS starts 
plunging to the BH at t ~ 90Af . Soon after this time, the 
NS is tidally disrupted; but the tidal disruption occurs 
near the ISCO and hence the material in the inner part is 
quickly swallowed into the BH (3rd panel). On the other 
hand, because of the outward angular momentum trans- 
fer, the material in the outer part of the NS forms a disk 
with the maximum density ~ 10 12 g/cm 3 (4th panel). 
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FIG. 2: (a) Evolution of baryon rest-mass located outside 
the apparent horizon for various grid settings. The plots for 
Ax — Afp /8 almost coincide and show that the results depend 
weakly on the values of L and £. On the other hand, the 
results depend systematically on Aa; (see text), (b) Evolution 
of averaged violation of Hamiltonian constraint for (iv)-(vii). 



However, mass of the disk is not large. Figure 2(a) 
shows evolution of baryon rest-mass located outside ap- 
parent horizons M r>rAH . We find that ~ 80% of the 
mass is swallowed into the BH in t ~ 130M ~ 2 ms, and 
swallowing continues after this time 33]. The value of 
M r>rAll depends systematically on Ax; we find that the 
results for (v)-(vii) at late times approximately obey a 
relation of convergence, i.e., M r>rAH (t) = a(t) + b(t)Ax n 
where a(t) and b(t) are functions of time. The order of 
convergence, denoted by n, is between 1st and 2nd or- 
der (i.e., 1 < n < 2). Least-square fitting gives a(t) at 
t = 180Af as w 0.19M* if we set n = 2 and as 0.28M, 
for n — 1 (see the solid circles in Fig. 2(a)). Thus, the 
true result should be between 0.19M* and 0.28M*. 

The adopted NS has corotating velocity field and fur- 
thermore its radius is larger than canonical values. In 
reality, the disk mass would be smaller than this value. 
Hence it is unlikely that a massive disk with ~ M Q is 
formed after merger of nonspinning BH of mass M > 
3M Q and canonical NS of mass w 1.4M Q and radius 
w 11-13 km, although a disk of mass of ~ 0.2-0.3Af Q 
may be formed for small BH mass and short GRBs of 
total energy ~ 10 49 ergs may be explained (e.g., [34j). 
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Figure 2(b) shows the evolution of averaged violation 
of the Hamiltonian constraint. For the average, rest-mass 
density is used as a weight (see [TEj]) and the integral is 
performed for the region outside apparent horizons. Fig. 
2(b) shows that the Hamiltonian constraint converges ap- 
proximately at 2nd order. This result is consistent with 
the fact that the region except for the vicinity of BH is 
followed with 2nd-order accuracy. 

To summarize, we have presented our first numerical 
results of fully GR simulation for merger of BH-NS bi- 
nary, focusing on the case that the BH is not spinning 
initially and the mass ratio q is fairly large as 0.4. It 
is found that even with such high value of q, the NS is 
tidally disrupted only for the orbit very close to ISCO 
and 80-90% of the mass element is quickly swallowed 
into the BH without forming massive disks. The results 
do not agree quantitatively with the prediction by the 
tidal approximation study. The reasons are: (1) In the 
tidal approximation, one describes NSs by the Newtonian 
gravity. In general relativity, gravity is stronger and the 
tidal disruption is less likely. (2) The time scale for angu- 



lar momentum transfer during tidal disruption near the 
ISCO is nearly as long as the plunging time scale deter- 
mined by gravitational radiation reaction and attractive 
force between two objects. Hence before the tidal disrup- 
tion completes, most of the material is swallowed. 

If the BH has a large spin, the final fate may be largely 
changed because of the presence of spin-orbit repulsive 
force. This force can weaken the attractive force between 
BH and NS and slow down the orbital velocity, resulting 
in smaller gravitational wave luminosity and longer ra- 
diation reaction time scale [35]. This effect may help 
massive disk formation. The study of spinning BH bina- 
ries is one of the next issues. The fate will also depend on 
EOS of NS [f| and mass of BH. Simulation with various 
EOSs and BH mass is also the next issue. 
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